library(sf)
library(dodgr)
library(geodist)
library(tidyverse)
library(sp)

# census data for 2020
AllLocalities = read.csv("Census2020.csv")
AllLocalitiesClean = subset(AllLocalities, AllLocalities$LONGITUD != "") ##Filtering out all
##non-locality data points (points without coordinates)

chd = "°"
chm = "'"
chs = "\""
AllLocalitiesClean$LONGITUD = ## converting degrees minutes seconds format to decimal format
  as.numeric(char2dms(AllLocalitiesClean$LONGITUD, chd = chd, chm = chm, chs = chs))
AllLocalitiesClean$LATITUD = ## converting degrees minutes seconds format to decimal format
  as.numeric(char2dms(AllLocalitiesClean$LATITUD, chd = chd, chm = chm, chs = chs))
AllLocalitiesClean$LONGITUD = round(AllLocalitiesClean$LONGITUD, 4)
AllLocalitiesClean$LATITUD = round(AllLocalitiesClean$LATITUD, 4) ## rounding to
## LAT and LONG to nearest 4th digit to be consistent with coop coordinates

LocalityShapeFile = st_as_sf(AllLocalitiesClean, coords = c(7,8), crs = 4326)

#census data for 2000
# chd = "°"
# chm = "'"
# chs = "\""
# AllLocalities2000 = read.csv("Census2000.csv")
# AllLocalitiesClean2000 = subset(AllLocalities2000, AllLocalities2000$long != "") ##Filtering out all
# ##non-locality data points (points without coordinates)
# AllLocalitiesClean2000 = subset(AllLocalitiesClean2000, AllLocalitiesClean2000$lat != "")
# AllLocalitiesClean2000 = subset(AllLocalitiesClean2000, select=c('entidad','nom_ent','mun', 'nom_mun', 'loc', 'nom_loc', 'long', 'lat', 'pobtot'))
# AllLocalitiesClean2000 = subset(AllLocalitiesClean2000, AllLocalitiesClean2000$long != ", LA (R")
# AllLocalitiesClean2000 = subset(AllLocalitiesClean2000, AllLocalitiesClean2000$long != ", LA")

# AllLocalitiesClean2000$long = ## converting degrees minutes seconds format to decimal format
#   as.numeric(char2dms(AllLocalitiesClean2000$long, chd = chd, chm = chm, chs = chs))
# AllLocalitiesClean2000$lat = ## converting degrees minutes seconds format to decimal format
#   as.numeric(char2dms(AllLocalitiesClean2000$lat, chd = chd, chm = chm, chs = chs))
# AllLocalitiesClean2000$long = round(AllLocalitiesClean2000$long, 4)
# AllLocalitiesClean2000$lat = round(AllLocalitiesClean2000$lat, 4)
# 
# LocalityShapeFile = st_as_sf(AllLocalitiesClean2000, coords = c(7,8), crs = 4326)
# LocalityShapeFile$POBTOT = LocalityShapeFile$pobtot

Coops = read.csv("CooperativePoints2022.csv") ## use file of cooperatives with location data here
## this file is a cleaned version of the raw CONAPESCA data. Much of this cleaning had to be
## done manually due to inconsistencies between census data and the RNPA. That is, we extract
## the longitude and latitude for each locality listed in the RNPA from the census. However, 
## many localities are not spelled the same way in the RNPA as they are in the census data. These 
## inconsistencies meant we often had to search the census manually to extract the proper coordinates 
## and attribute them to the localities listed in the RNPA. We were not able to find some localities
## listed in the RNPA in the census. These are included in the number of cooperatives that have
## missing location data. 
Coops = subset(Coops, Coops$MOTIVO != "DUPLICIDAD") ## get rid of duplicates
Coops = subset(Coops, Coops$LOCALIDAD != "ISLA MUJERES") ## get rid of coops on islands that are not connected to cities via roads
Coops = subset(Coops, Coops$LOCALIDAD != "COZUMEL")
Coops = subset(Coops, Coops$UNIDAD.ECONÓMICA != "X") ## get rid of coops without names
Coops = subset(Coops, Coops$UNIDAD.ECONÓMICA != "")
Coops = subset(Coops, Coops$LAT != "") # get rid of cooperatives that do not have location data


CoopsShapeFile = st_as_sf(Coops, coords = c("long", "lat"), crs = 4326) ## Converting coop data to 

## shape file 
##Coops = subset(CoopsClean, CoopsClean$MOTIVO != "DUPLICIDAD")


# change working directory to the roads folder

mxroads = st_read("red_vial.shp")

mxroads = subset(mxroads, mxroads$TIPO_VIAL != "Andador")
mxroads = subset(mxroads, mxroads$TIPO_VIAL != "Peatonal")
mxroads = subset(mxroads, mxroads$VELOCIDAD != "N/A")
mxroads = subset(mxroads, mxroads$VELOCIDAD != "5")

mxroads$ID = 1:nrow(mxroads)


## weighting road distance according to speed limit so output is travel time
colnm <- "VELOCIDAD"
wts <- data.frame (
  name = "custom",
  way = unique(mxroads[[colnm]]),
  value = c (80000, 60000, 50000, 40000, 20000, 30000, 70000, 10000, 90000, 110000, 100000)
)

#establishing our weighted network
g = weight_streetnet(mxroads, wt_profile = wts, type_col = "VELOCIDAD", id_col = "ID")
g$d = g$d_weighted

# getting subsets of census localities that we will measure distance to
TenK = filter(LocalityShapeFile, POBTOT >= 10000)
TwentyK = filter(LocalityShapeFile, POBTOT >= 20000)
ThirtyK = filter(LocalityShapeFile, POBTOT >= 30000)
FortyK = filter(LocalityShapeFile, POBTOT >= 40000)
FiftyK = filter(LocalityShapeFile, POBTOT >= 50000)
SixtyK = filter(LocalityShapeFile, POBTOT >= 60000)
SeventyK = filter(LocalityShapeFile, POBTOT >= 70000)
EightyK = filter(LocalityShapeFile, POBTOT >= 80000)
NinetyK = filter(LocalityShapeFile, POBTOT >= 90000)
OneHundredK = filter(LocalityShapeFile, POBTOT >= 100000)
OneTenK = filter(LocalityShapeFile, POBTOT >= 110000)
OneTwentyK =  filter(LocalityShapeFile, POBTOT >= 120000)
OneThirtyK = filter(LocalityShapeFile, POBTOT >= 130000)
OneFortyK = filter(LocalityShapeFile, POBTOT >= 140000)
OneFiftyK = filter(LocalityShapeFile, POBTOT >= 150000)
OneSixtyK = filter(LocalityShapeFile, POBTOT >= 160000)
OneSeventyK = filter(LocalityShapeFile, POBTOT >= 170000)
OneEightyK = filter(LocalityShapeFile, POBTOT >= 180000)
OneNinetyK = filter(LocalityShapeFile, POBTOT >= 190000)
TwoHundredK = filter(LocalityShapeFile, POBTOT >= 200000)
TwoFiftyK = filter(LocalityShapeFile, POBTOT >= 250000)
ThreeHundredK = filter(LocalityShapeFile, POBTOT >= 300000)
FourHundredK = filter(LocalityShapeFile, POBTOT >= 400000)
FiveHundredK = filter(LocalityShapeFile, POBTOT >= 500000)
Million= filter(LocalityShapeFile, POBTOT >= 1000000)

#10k
from = st_coordinates(CoopsShapeFile) #points we are measuring from (cooperatives)
to10k = st_coordinates(TenK) #points we are measuring to (localities of at least 10,000)
d10k = dodgr_dists(graph = g, from = from, to = to10k)
df10k = apply(d10k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df10k) #adding travel time to nearest locality of at least 10,000 to our cooperative dataset 
names(Coops)[23] = "time10k"

#20k
to20k = st_coordinates(TwentyK)
d20k = dodgr_dists(graph = g, from = from, to = to20k)
df20k = apply(d20k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df20k)
names(Coops)[24] = "time20k"

#30k
to30k = st_coordinates(ThirtyK)
d30k = dodgr_dists(graph = g, from = from, to = to30k)
df30k = apply(d30k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df30k)
names(Coops)[25] = "time30k"

#40k
to40k = st_coordinates(FortyK)
d40k = dodgr_dists(graph = g, from = from, to = to40k)
df40k = apply(d40k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df40k)
names(Coops)[26] = "time40k"

# #50k
to50k = st_coordinates(FiftyK)
d50k = dodgr_dists(graph = g, from = from, to = to50k)
df50k = apply(d50k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df50k)
names(Coops)[27] = "time50k"


# 50k with name of nearby population
# for (i in 1:5427) {
#   coop = matrix(1, nrow = 232, ncol = 5)
#   coop = data.frame(coop)
#   from = st_coordinates(CoopsShapeFile[i,])
#   for(j in 1:232) {
#     to = st_coordinates(FiftyK[j,])
#     coop$X1[j] = dodgr_dists(graph = g, from = from, to = to)
#     coop$X2[j] = FiftyK$NOM_LOC[j]
#     coop$X3[j] = FiftyK$NOM_MUN[j]
#     coop$X4[j] = FiftyK$NOM_ENT[j]
#     coop$X5[j] = FiftyK$POBTOT[j]
#   }
#   CoopsShapeFile$Distance50K[i] = min(as.numeric(coop$X1))
#   for (k in 232) {
#     if (coop$X1[k] == min(coop$X1)) {
#       CoopsShapeFile$L50K[i] = FiftyK$NOM_LOC[k]
#       CoopsShapeFile$M50K[i] = FiftyK$NOM_MUN[k]
#       CoopsShapeFile$E50K[i] = FiftyK$NOM_ENT[k]
#       CoopsShapeFile$P50K[i] = FiftyK$POBTOT[k]
#       break
#     }
#   }
# }

# #60k
to60k = st_coordinates(SixtyK)
d60k = dodgr_dists(graph = g, from = from, to = to60k)
df60k = apply(d60k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df60k)
names(Coops)[28] = "time60k"

# #70k
to70k = st_coordinates(SeventyK)
d70k = dodgr_dists(graph = g, from = from, to = to70k)
df70k = apply(d70k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df70k)
names(Coops)[29] = "time70k"

# #80k
to80k = st_coordinates(EightyK)
d80k = dodgr_dists(graph = g, from = from, to = to80k)
df80k = apply(d80k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df80k)
names(Coops)[30] = "time80k"
# 
# #90k
to90k = st_coordinates(NinetyK)
d90k = dodgr_dists(graph = g, from = from, to = to90k)
df90k = apply(d90k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df90k)
names(Coops)[31] = "time90k"

#100k
to100k = st_coordinates(OneHundredK)
d100k = dodgr_dists(graph = g, from = from, to = to100k)
df100k = apply(d100k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df100k)
names(Coops)[32] = "time100k"

# #110k
to110k = st_coordinates(OneTenK)
d110k = dodgr_dists(graph = g, from = from, to = to110k)
df110k = apply(d110k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df110k)
names(Coops)[33] = "time110k"
# 
# #120k
to120k = st_coordinates(OneTwentyK)
d120k = dodgr_dists(graph = g, from = from, to = to120k)
df120k = apply(d120k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df120k)
names(Coops)[34] = "time120k"

# 
# for (i in 1:5427) {
#   coop = matrix(1, nrow = 129, ncol = 5)
#   coop = data.frame(coop)
#   from = st_coordinates(CoopsShapeFile[i,])
#   for(j in 1:129) {
#     to = st_coordinates(OneTwentyK[j,])
#     coop$X1[j] = dodgr_dists(graph = g, from = from, to = to)
#     coop$X2[j] = OneTwentyK$NOM_LOC[j]
#     coop$X3[j] = OneTwentyK$NOM_MUN[j]
#     coop$X4[j] = OneTwentyK$NOM_ENT[j]
#     coop$X5[j] = OneTwentyK$POBTOT[j]
#   }
#   CoopsShapeFile$Distance120K[i] = min(as.numeric(coop$X1))
#   for (k in 129) {
#     if (coop$X1[k] == min(coop$X1)) {
#       CoopsShapeFile$L120K[i] = OneTwentyK$NOM_LOC[k]
#       CoopsShapeFile$M120K[i] = OneTwentyK$NOM_MUN[k]
#       CoopsShapeFile$E120K[i] = OneTwentyK$NOM_ENT[k]
#       CoopsShapeFile$P120K[i] = OneTwentyK$POBTOT[k]
#     }
#   }
# }

# 
# #130k
to130k = st_coordinates(OneThirtyK)
d130k = dodgr_dists(graph = g, from = from, to = to130k)
df130k = apply(d130k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df130k)
names(Coops)[35] = "time130k"
# 
# #140k
to140k = st_coordinates(OneFortyK)
d140k = dodgr_dists(graph = g, from = from, to = to140k)
df140k = apply(d140k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df140k)
names(Coops)[36] = "time140k"

#150k
to150k = st_coordinates(OneFiftyK)
d150k = dodgr_dists(graph = g, from = from, to = to150k)
df150k = apply(d150k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df150k)
names(Coops)[37] = "time150k"

#160k
to160k = st_coordinates(OneSixtyK)
d160k = dodgr_dists(graph = g, from = from, to = to160k)
df160k = apply(d160k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df160k)
names(Coops)[38] = "time160k"

#170k
to170k = st_coordinates(OneSeventyK)
d170k = dodgr_dists(graph = g, from = from, to = to170k)
df170k = apply(d170k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df170k)
names(Coops)[39] = "time170k"

#180k
to180k = st_coordinates(OneEightyK)
d180k = dodgr_dists(graph = g, from = from, to = to180k)
df180k = apply(d180k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df180k)
names(Coops)[40] = "time180k"

#190k
to190k = st_coordinates(OneNinetyK)
d190k = dodgr_dists(graph = g, from = from, to = to190k)
df190k = apply(d190k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df190k)
names(Coops)[41] = "time190k"

#200k
to200k = st_coordinates(TwoHundredK)
d200k = dodgr_dists(graph = g, from = from, to = to200k)
df200k = apply(d200k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df200k)
names(Coops)[42] = "time200k"

#250k
to250k = st_coordinates(TwoFiftyK)
d250k = dodgr_dists(graph = g, from = from, to = to250k)
df250k = apply(d250k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df250k)
names(Coops)[43] = "time250k"

#300k
to300k = st_coordinates(ThreeHundredK)
d300k = dodgr_dists(graph = g, from = from, to = to300k)
df300k = apply(d300k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df300k)
names(Coops)[44] = "time300k"

#400k
to400k = st_coordinates(FourHundredK)
d400k = dodgr_dists(graph = g, from = from, to = to400k)
df400k = apply(d400k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df400k)
names(Coops)[45] = "time400k"

#500k
to500k = st_coordinates(FiveHundredK)
d500k = dodgr_dists(graph = g, from = from, to = to500k)
df500k = apply(d500k, 1, FUN=min, na.rm=TRUE) %>% as.data.frame()
Coops = cbind(Coops, df500k)
names(Coops)[46] = "time500k"



